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Abstract 

We  consider  the  problem  of  tracing  algebraic  curves  by  computer, 
using  a  numerical  technique  augmented  by  symbolic  computations.  In 
particular,  all  singularities  are  analyzed  correctly.  The  methods  pre¬ 
sented  find  application  in  solid  modeling  and  robotics. 


1  Introduction 

In  this  paper  we  discuss  preliminary  results  for  tracing  algebraic  curves. 
Planar  algebraic  curves  of  the  form  /(*,  y)  =  0  are  considered,  as  are  space 
curves  that  are  the  intersection  of  two  algebraic  surfaces,  f(x,y,z)  =  0  and 
g(x,y,z)  —  0.  The  scenario  is  as  follows:  We  are  given  a  point  p  on  the  curve, 
and  a  direction  of  traversal.  We  wish  to  trace  succeeding  curve  points  on 
the  same  branch,  and  we  would  like  to  trace  them  through  singularities. 

A  reliable  solution  to  this  problem  has  immediate  applications  to  solid 
modeling  and  geometric  design.  For  example,  the  well-known  Boolean  oper¬ 
ations  on  solids  require  tracing  surface  intersections,  for  the  purpose  of  de¬ 
termining  the  surface  of  the  resulting  solid  [9).  Here,  algebraic  space  curves 
arise  when  the  intersecting  solids  are  bounded  by  algebraic  faces.  If  the 
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surface  is  composed  of  parametric  patches,  e.g.,  rational  B-spline  surfaces, 
then  planar  algebraic  curves  can  be  obtained  [6]. 

Since  an  extensive  amount  of  curve  tracing  is  required  for  each  modeling 
operation,  it  is  advisable  to  pay  attention  to  efficiency  as  well  as  reliability. 
For  this  reason,  precise  methods  such  as  the  cylindrical  algebraic  decompo¬ 
sition  due  to  Collins  and  its  variants  [4,11]  have  not  been  considered  here. 
This  does  not  imply  that  these  presently  very  compute  intensive  methods 
must  remain  of  theoretical  interest  c  nly,  but  significantly  more  work  is  re¬ 
quired  before  we  can  accurately  gauge  whether  in  the  specialized  context  of 
solid  modeling  versions  of  these  algorithms  exist  that  can  be  used  without 
serious  efficiency  degradations. 

As  point  of  departure,  we  use  a  straightforward  numerical  method  that 
approximates  the  curve  locally  by  its  truncated  Taylor  series,  and  then  per¬ 
forms  a  Newton  iteration  to  correct  the  accumulated  error.  This  approxi¬ 
mation  has  many  interesting  properties.  For  instance,  the  Taylor  series  is 
a  special  case  of  the  Puiseux  series  that  can  be  used  to  approximate  the 
curve  at  singular  points  as  well  as  at  regular  points.  Thus,  there  is  a  basis 
for  developing  a  uniform  framework  for  studying  approximations  to  alge¬ 
braic  curves.  Note,  however,  that  there  are  alternative  curve  approximations 
based  on  special  classes  of  polynomials  that  offer  different  advantages,  [7], 
and  more  study  is  needed  before  the  relative  merits  can  be  fully  appreciated. 

In  many  cases  the  numerical  procedure  suffices  and  copes  acceptably  well 
with  certain  singularities,  e.g.,  with  normal  crossings  and  with  tacnodes  that 
are  not  very  complicated.  It  is  not  fully  reliable,  however,  and  will  fail  at 
cuspidal  singularities.  For  this  reason  it  is  augmented  by  a  mapping  tech¬ 
nique  that  exploits  the  fact  that  by  a  suitable  birational  map  any  singularity 
can  be  resolved.  That  is,  such  a  map  will  transform  the  singular  point  into 
one  or  more  nonsingular  ones,  while  not  creating  new  singularities.  Fortu¬ 
nately,  suitable  maps  can  be  found  easily  in  the  planar  case.  In  the  space 
curve  case  the  situation  is  not  so  simple,  and  more  work  is  required  to  find 
attractive  algorithms. 

Even  in  the  planar  case  a  number  of  details  must  be  addressed  before 
curve  desingularization  can  be  automated.  These  include  finding  reliably  the 
locus  of  the  singularity,  controlling  numerical  inaccuracies  that  arise  from  the 
various  desingularization  maps,  and  establishing  the  correct  correspondence 
of  orientation  between  the  curve  and  its  transformations. 

We  structure  this  paper  as  follows:  After  explaining  concepts  and  nota- 
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tion  in  Section  2,  we  devote  Section  3  to  a  description  of  the  simple  numerical 
procedure  for  following  curves.  In  Section  4,  we  explain  the  correspondence 
between  the  Taylor  series  and  the  notion  of  a  place  of  a  curve,  which  yields  a 
straightforward  method  for  extending  the  numerical  procedure  so  as  to  cope 
with  low  order  singularities  consistently.  The  extension  is  of  limited  value, 
however,  as  it  involves  solving  systems  of  polynomial  equations  whose  degree 
depends  on  the  order  of  the  singularity  analyzed.  Section  5  concentrates  on 
desingularizing  planar  algebraic  curves  and  discusses  how  problems  such  as 
orientation  correspondence  can  be  solved.  Examples  of  various  types  of  sin¬ 
gularities  are  also  given.  Section  6  then  discusses  the  desingularization  of 
space  curves.  Sections  3,  4  and  5  summarize  work  reported  in  [2,8], 


2  Concepts  and  Notation 

We  consider  algebraic  space  curves  given  as  the  intersection  of  two  algebraic 
surfaces  /(x,y,z)  =  0  and  g[x,y,z)  —  0,  where  /  and  g  are  polynomials  in 
x,  y,  and  z.  This  is  not  the  most  general  definition  of  space  curves;  certain 
curves  require  intersecting  more  than  two  surfaces  in  order  to  exclude  extra¬ 
neous  components.  The  partial  derivatives  of  /  are  written  by  subscripting, 
e.g.,  /*  =  df  /dx,  f xy  =  d2f/[dxd y),  and  so  on.  Since  /  and  g  are  analytic, 

fxy  ~  fyx  etc. 

Vectors  and  vector  functions  are  denoted  by  bold  letters.  The  inner 
product  of  two  vectors  a  and  b  is  the  scalar  a  •  b.  The  length  of  the  vector 
a  is  |a|  =  y/a  •  a.  The  cross  product  of  the  vectors  is  the  vector  a  x  b. 

The  gradient  of  the  surface  /  at  the  point  p  =  (x,y,z)  is  the  vector 
V/  =  (/u/y,/*),  where  the  partials  are  evaluated  at  p.  If  not  zero,  it  is  a 
vector  that  is  normal  to  the  surface  at  p.  A  point  p  =  (x,y,  z)  is  regular  on 
/  if  the  gradient  of  /  at  p  is  not  null;  otherwise  the  point  is  singular.  The 
Hessian  of  /  at  point  p  is  the  tensor 

(fxx  fxy  fxz  \ 
fyx  fyy  fyx  I 
f xx  fxy  fxx  J 

where  the  partials  are  evaluated  at  p. 

The  intersection  curve  of  /  and  g  is  denoted  r(s)  and  is  considered  a 
vector  function  of  the  argument  s.  In  Section  3,  we  will  determine  an  ap¬ 
proximation  of  r  in  which  s  is  the  arc  length  measured  from  some  initial  point 
on  the  curve.  As  usual,  the  derivatives  of  r(s)  are  denoted  r',  r" . . 
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A  point  p  of  the  intersection  curve  r(e)  is  regular  if  p  is  regular  on  both 
/  and  g  and  if  the  gradients  V/  and  Vg  are  linearly  independent.  That  is, 
the  surfaces  are  not  singular  at  p  and  intersect  transversally.  A  point  p  is 
singular  on  r(s)  for  one  of  the  following  reasons: 

1.  The  gradients  V /  and  Vg  are  nonzero  and  linearly  dependent. 

2.  One  of  the  gradients,  say  Vg  is  zero,  but  the  other  is  not. 

3.  The  gradients  V/  and  Vg  are  both  zero. 


We  note  that  Cases  1  and  2  do  not  differ  in  substance. 

The  initial  form  of  /  is  the  polynomial  formed  by  all  terms  of  lowest 
order  in  /.  For  example,  the  initial  form  of  z2  -  2a:  +  y2  +  z1  is  —  2s.  The 
initial  form  approximates  the  surface  in  the  neighborhood  of  the  origin.  In 
the  above  example,  the  surface  “looks  like”  the  plane  z  =  0  near  the  origin. 
This  plane  is  the  tangent  plane  to  /  at  the  origin.  In  particular,  if  p  is  the 
origin,  then  p  is  regular  on  /  if  the  initial  form  of  /  is  linear.  Otherwise  p 
is  singular. 

At  each  point  p  =  r(s)  of  a  space  curve,  an  intrinsic  coordinate  system  is 
provided  by  the  orthonormal  triad.  The  triad  consists  of  three  perpendicular 
unit  vectors,  namely  the  tangent  t,  the  principal  normal  h,  and  the  binormal 
b.  The  Frenet-Serret  formulas  relate  the  triad  to  arc  length,  curvature,  and 
torsion  of  the  curve  at  p.  With  s  the  arc  length,  p  the  radius  of  curvature, 
and  r  the  radius  of  torsion  we  have 


dt 

ds 


db  1,  dh  1,  1 

—  =  --h,  _  =  -b--t 

ds  r  ds  r  p 


A  planar  algebraic  curve  is  given  by  its  implicit  equation  /(z,y)  =  0, 
/  a  polynomial  in  z  and  y.  Equivalently,  we  can  think  of  the  curve  as  the 
intersection  of  the  surfaces  /  and  z  =  0.  Note  that  /  is  then  a  cylinder  with 
the  base  line  /(z,  y)  =  0.  All  concepts  explained  above  can  therefore  be 
transferred  to  the  planar  case. 


3  Numerical  Tracing 

For  space  curves,  the  simplest  situation  arises  when  tracing  the  curve  in 
a  neighborhood  in  which  both  surfaces  are  nonsingular  and  intersect  each 
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other  transversally.  This  means  that  the  gradients  of  /  and  g  do  not  van¬ 
ish  along  the  curve  and  are  linearly  independent  vectors.  In  this  case  we 
formulate  a  system  of  linear  equations  from  which  to  obtain  the  local  ap¬ 
proximation  to  the  intersection  curve  at  the  point  p. 

3.1  The  Basic  Method  for  Space  Curves 

Since  the  curve  r(s)  must  satisfy  both  /  and  g  identically,  each  coefficient 
in  the  Taylor  series  of  /(r(s))  and  g(r(s))  will  be  zero.  Let  p  =  r(0)  be  a 
point  on  the  intersection.  Then 

/(*(«))  =  /( 0)  +  «V/  ■  r'(0)  +  j [V/  •  r"(0)  +  r'(0)  •  Hf  •  r'(0)]  +  •  •  ■ 
and  similarly, 

g(r(s))  =  g(0)  +  sVg  •  r'(0)  +  y  [VS  •  r"(0)  +  r'(0)  •  Ha  ■  r'(0)]  +  •  •  • 
This  leads  to  the  following  system  of  equations,  for  m  =  1,2, . . .: 

V/(p).r(”l(0)  B,,m 

V9(p)tW(0)  =  (1) 

The  quantities  Bi<m  and  B2,m  are  expressions  in  the  partial  derivatives  of  / 
and  g  and  the  lower-order  derivatives  of r.  For  m=lwe  have 

B-  i  =  $2,1  =  0 

and  for  m  =  2 

$1,2  =  -r'-Hf-r' 

$2,2  =  -r'  •  Hg  •  r' 

For  higher  values  of  m  the  expressions  $,,m  are  more  complex. 

With  the  assumption  of  independent  gradients  at  r(0),  the  solution  to 
the  system  has  the  form 

<*mV/  +  0mVg  +  TfmV/  X  V? 

7m  can  be  chosen  arbitrarily,  and  the  other  coefficients  satisfy  the  nonsin¬ 
gular  system 

/V/.V/  V/-VjWam\_/  $i,m\  (2) 

\Vg-Vf  Vg-VgJ\/3mJ 
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System  (l)  is  underdetermined.  When  using  it  to  approximate  the  in¬ 
tersection  curve,  the  choices  of  7 m  relate  to  the  ability  to  obtain  equivalent 
approximations,  e.g.,  parameterized  by  cs  instead  of  a,  where  e  ^  0  is  a  con¬ 
stant.  By  using  the  Frenet-Serret  formulas,  we  can  interpret  these  choices. 
We  set  r'  =  dr/ds  =  t,  to  obtain 


For  m  =  1  we  have  Biti  =  £2,1  =  0,  hence  oc\  =  =  0.  We  let 

r,_  V/xVg 
IV/  x  Vff| 

be  the  solution  of  the  system,  so  that  the  parameter  s  corresponds  to  the 
arc  length  of  the  intersection  curve. 

For  m  =  2  we  choose  72  =  0-  This  implies  that  r"  is  orthogonal  to  r', 
and  that  from  the  system  solution  r"  =  ct^V f  +  fcVg  =  h/p  both  the 
principal  normal  h  and  the  radius  of  curvature  p  is  determined. 

For  m  =  3  we  get  73  =  —  (r"  •  r"),  since  both  b  and  h  are  orthogonal  to 
t.  This  determines  both  r'"  and  r. 


The  curve  approximation  so  determined  can  be  used  in  a  neighborhood 
of  r(0)  whose  size  can  be  estimated  by  the  magnitude  of  the  second  and 
third  order  terms.  If  the  ratio  of  $s|rm|/6  to  |r  +  Sr'  +  62t"/2\  is  small,  then 
the  higher  order  derivatives  contribute  little  and  we  have  not  deviated  too 
much  from  the  true  intersection.  By  halving  or  doubling  a  standard  distance 
repeatedly,  the  stepping  size  can  be  adjusted  according  to  curvature  and 
torsion.  It  is  necessary  to  establish  a  minimum  stepping  size  since  near¬ 
singular  curves  can  have  areas  of  arbitrarily  high  curvature  where  repeated 
halving  might  lead  to  unacceptable  running  times. 

At  the  end  of  the  current  approximation  to  r  a  point  Po  is  reached  that 
is  near  the  curve  of  intersection  but  not  on  it.  Beginning  with  this  point, 
we  find  a  sequence  of  points  Pi,P2,...  that  converges  to  a  point  p  on  the 
curve,  using  Newton’s  method.  With  Pk+i  =  P*  +  A*,  we  want  to  solve  the 
system 


V/(P*)-A*  =  -/(P*) 

Vp(P*)  •  A*  =  -g(Pk)  (3) 
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to  obtain  A*  rs  sP'*.  Assuming  linearly  independent  gradients, 

Pk  =  <**V/  +  Pk^g  +  7*t, 

where  t  is  orthogonal  to  both  gradients  evaluated  at  P*.  Since  a  change 
in  the  direction  of  t  does  not  change  the  values  of  /  and  g  appreciably, 
we  set  7*  =  0,  thereby  obtaining  a  unique  solution  for  A*.  We  then  set 
Pfc+i  =  Pk  +  A*. 

After  the  point  p  is  found  with  acceptable  accuracy,  a  new  approximation 
of  r(s)  centered  at  p  is  determined. 

3.2  The  Planar  Case 

The  planar  curve  /  (*,  y)  =  0  arises  as  intersection  of  the  cylinder  /  (*,  y)  =  0 
and  the  plane  z  =  0.  It  can  be  traced  this  way  as  a  space  curve.  In  the 
system  (1)  the  second  equation  specializes  to  r =  0.  Here  denotes 
the  z  component  of  the  mth  derivative  of  r.  Moreover,  since  all  partial 
derivatives  by  z  of  /  are  zero,  the  first  equation  takes  the  form 

/*rim)  +  /vr^)  =  Cm 

Thus,  there  is  no  difference  between  considering  the  intersection  curve  r  or 
the  planar  curve. 

As  before,  choosing  7  =  0  for  the  Newton  iteration  means  that  we  ap¬ 
proach  the  curve  along  the  local  normal  direction.  An  implementation  could 
be  specialized,  but  there  appears  to  be  no  significant  penalty  for  tracing  the 
curve  in  space. 

3.3  Implementation 

The  numerical  tracing  procedure  has  been  implemented  in  Fortran  by  R. 
Lynch  on  a  VAX  8600.  With  minor  modifications,  it  has  then  been  ported 
to  a  Symbolics  Lisp  machine.  Figures  3.1  through  3.4  show  some  examples 
of  curve  traces  so  obtained.  The  planar  curves  have  been  traced  as  the 
intersection  of  f(x,  y)  =  0  with  z  =  0. 

In  our  experience,  nodal  singularities  cause  no  problems  as  long  as  the 
tangent  directions  of  the  intersecting  branches  are  sufficiently  separated. 
Since  the  curve  orientation  may  reverse  at  singularities  (c.f.  Subsection  5.4 
below),  the  tracing  program  must  be  augmented  so  as  to  maintain  consistent 
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tangent  direction.  However,  the  program  cannot  trace  through  cuspidal 
singularities.  Many  tacnodes  are  handled  reliably,  but  inflections  at  the 
singularity  are  not  recognized.  Thus  both  the  curve  f\  =  y2  -  x4  -  y4  =  0, 
shown  in  Figure  3.3,  and  the  curve  ft  =  y2  -  z6  -  y6  =  0,  shown  in  Figure 
3.4,  are  traced  as  if  they  have  two  real  components  tangentially  meeting  at 
the  origin.  While  this  is  correct  for  / 1,  it  is  not  correct  for  fi  which  consists 
of  a  single  real  component  with  the  two  branches  at  the  origin  each  having 
a  point  of  inflection. 


4  Algebraic  Extensions  of  the  Method 

We  derived  the  equations  (1)  based  on  the  assumption  that  the  two  surfaces 
intersect  transversally  and  are  not  singular  at  the  point  r(0)  of  interest. 
Clearly,  the  radius  of  convergence  of  the  power  series  about  such  a  point 
cannot  include  any  singular  points  of  the  curve.  Nonetheless,  the  system  of 
equations  remains  valid  even  when  we  are  at  a  singular  curve  point.  The 
reason  for  this  is  that  Taylor’s  theorem  is  a  special  case  of  more  general 
theorems. 

Informally,  a  place  of  the  planar  curve  /( x,y)  —  0  is  a  pair  of  power 
series 

=  y{s)  =  ^2bksk 

k> 0  *>0 

such  that  f(x[s),y(s))  is  identically  zero.  The  center  of  a  place  is  the  point 
(z(0),y(0))  on  the  curve.  Newton’s  Theorem  states  that  centered  at  every 
curve  point  there  is  at  least  one  place  of  /.  Likewise,  we  define  a  place  of 
the  space  curve  r  as  the  triple  (i(s),y(s), z(s)).  Since  every  space  curve  is 
the  birational  image  of  a  planar  curve  [13,14],  there  is  at  least  one  place 
centered  at  every  point  of  the  space  curve. 

The  connection  between  the  notion  of  place  and  the  equation  system  (l) 
in  the  previous  section  is  established  as  follows.  Centered  at  the  point  p  we 
consider  the  place 

r(s)  =  ^2{ak,h,Ck)sk 

k>  0 

where  p  =  r(0).  The  derivative  of  the  place  is  defined  by 

As)  =  ^(ak,bk,Ck)ksk-k 
t>  l 
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Higher  order  derivatives  are  defined  analogously. 

Since  p  is  assumed  to  be  on  the  intersection  curve  of  /  and  g ,  we  know 
that  /(r(s))  and  g(r(s))  are  identically  zero,  from  which  a  system  of  equa¬ 
tions  is  obtained  for  m  =  1, 2, 3, . . . 

Kx<m  =  0 

Ki.r,  =  o  (4) 

Here  is  the  coefficient  of  sm  in  the  power  series  /(r(s))  and  Jf2>m  is 
the  coefficient  of  sm  in  the  power  series  p(r(s)).  This  leads  to  the  following 

Theorem  For  all  m  >  1  the  equation  V/r^m^(0)  =  B iiTO  is  equivalent 
to  the  equation  m\Ki<m  =  0. 

The  analogous  statement  holds  for  g  and  A^.m-  The  proof  is  by  induction 
on  the  terms  of  / ;  see  [8]  for  details.  How  the  series  are  obtained  from  system 
(4)  is  illustrated  by  an  example. 

Consider  the  cylinders  /  =  x2  +  y2  +  2x  =  0  and  g  =  i2  +  z2  +  4x  =  0. 
Their  intersection  is  an  irreducible  space  curve  of  degree  4  with  a  singular 
point  at  the  origin.  At  the  singular  point  we  have  the  following  equations: 

a\  =  0 

ai  =  0 

dj  +  2a2  +  b\  =  0 

cj  +  4a2  +  oj  =  0 

2aia2  +  2a3  +  26162  =  0 

2djd2 -(- 4d3  +  2cjc2  —  0 

d2  +  2ajd3  -f-  2d4  +  62  +  26163  =  0 

o2  +  2did3  +  4d4  +  c\  +  2ciC3  =  0 


One  of  the  solutions  to  this  system  is 


x(s) 

y(s) 

z(s) 


V2 s - ^=s8  ■ 

2v/2 


2s  -  -s3 
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In  principle,  this  approach  can  be  used  to  extend  the  tracing  method  of 
Section  3  so  as  to  handle  singular  points,  but  it  may  become  computationally 
expensive.  Efficient  strategies  for  solving  these  equations  may  exist.  For 
example,  one  can  always  choose  the  coefficients  of  one  series,  say  for  x(s), 
such  that  |a*|  =  1  for  one  specific  k,  and  all  other  coefficients  are  zero  [14]. 


5  Planar  Algebraic  Curves 

The  problem  of  tracing  a  curve  reliably  through  any  singularity  is  partially 
solved  by  the  following  theorem  from  algebraic  geometry,  e.g.  [1]: 

Theorem  Any  given  algebraic  plane  curve  can  be  transformed,  by  a 
birational  transformation,  into  a  curve  devoid  of  singularities. 

The  proof  proceeds  by  an  inductive  argument  that  builds  up  the  required 
birational  transformation  through  a  sequence  of  elementary,  quadratic  trans¬ 
formations.  It  is  easy  to  understand  that  these  transformations  resolve  or¬ 
dinary  singularities,  but  how  progress  is  made  on  irregular  singularities  is 
more  subtle,  and  we  will  not  discuss  it  here.  Different  proofs  of  the  theorem 
are  found  in,  e.g.,  [1,10,13,14].  We  restrict  our  attention  to  those  transfor¬ 
mation  properties  that  are  needed  in  order  to  understand  how  to  derive  an 
algorithm  from  the  theorem. 

5.1  Desingularization 

We  map  a  planar  curve  /(*,  y)  =  0  to  a  curve  g( x\,yi) 
transformation 

y 

xi  =  x  yi=  - 
x 

If  /  has  a  singular  point  of  order  m  at  the  origin,  then 

f{*i,x  iyi)  =  *rs(*i.yi) 

We  call  f(xi,xiyi)  the  total  transform,  and  y(xi,yi)  the  proper  transform 
of  /.  The  yi-axis  is  called  the  exceptional  line.  Figures  5.1  and  5.2  show 
two  examples  of  a  curve  /  and  its  proper  transform. 

Intuitively  speaking,  applying  the  quadratic  transformation  separates 
intersecting  curve  branches  that  have  different  tangent  directions.  To  ap¬ 
preciate  this,  note  that  the  line  y  —  mx  =  0  is  transformed  to  the  line 


=  0  by  the  quadratic 

(5) 
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j/i  -  m  =  0.  Moreover,  all  points  (x,  y )  of  the  x-y  plane  with  r  ^  0  are  in 
1  —  1  correspondence  with  points  ( x,y/x )  of  the  zi-yi  plane.  A  point  (0,y) 
of  the  x-y  plane  with  y  ^  0  is  mapped  to  infinity  in  the  zi-yi  plane,  and  the 
origin  of  the  x-y  is  mapped  to  the  exceptional  line.  The  effect  is  that  the 
singular  point  is  “blown  up”  to  the  line  xi  =  0,  and  that  the  branches  of  / 
at  the  origin  are  separated  or,  in  a  precise  sense,  made  less  singular.  The 
proof  of  the  theorem  shows  that  after  a  finite  number  of  quadratic  transfor¬ 
mations  all  singularities  are  removed.  We  wish  to  trace  through  a  singular 
point  as  follows: 

1.  We  trace  the  curve  /  using  the  basic  method  of  Section  3. 

2.  When  approaching  a  singularity,  we  notice  at  some  point  p  that  the 
determinant  of  the  system  becomes  too  small.  We  then  locate  the  sin¬ 
gular  point  as  described  below,  and  move  it  to  the  origin  by  translating 
the  coordinate  system. 

3.  Now  the  quadratic  transformation  is  applied  yielding  the  proper  trans¬ 
form  g. 

4.  We  traverse  g  beginning  at  the  point  pi  corresponding  to  p,  until  we  are 
past  the  singularity  of  /  and  the  system  determinant  is  large  enough 
to  continue  traversing  /  accurately. 

Note  that  we  may  have  to  traverse  recursively  iterated  transforms  of  /, 
since  the  applied  quadratic  transform  may  not  have  fully  desingularized  the 
corresponding  branch  of  g.  Moreover,  care  must  be  exercised  in  correlating 
the  orientation  of  g  and  of  /  to  maintain  proper  traversal  direction. 

5.2  Locating  the  Singularity 

When  approaching  a  singular  point  p,  the  partial  derivatives  fz  and  /„  of 
/  vanish.  If  the  singularity  has  higher  order,  then  higher  order  partial 
derivatives  also  vanish. 

The  singular  point  is  defined  as  the  intersection  of  the  curves  /  =  0, 
fz  =  0,  and  fv  =  0.  When  traversing  the  curve  /,  we  have  approached 
the  singularity  to  a  point  Pq  for  which  the  partial  fx  and  fv  drop  in  value 
below  a  threshold  p.  We  use  a  Newton  iteration  to  construct  a  sequence  of 
point  approximations  P,  that  converges  to  the  singularity.  The  iteration  is 
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governed  by  the  following  system: 


(  U(Pi) 

fv(Pi)  ' 

1/f  \ 

/zz(P.) 

hv(Pi) 

f  ) 

^  /zy(Pi) 

fvv(P')  j 

l  \°v  J 

m )  ^ 
/,(«) 
/.(«)  j 


whose  solution  determines  the  next  approximation  to  the  intersection  as 
Pt+i  =  Pi  +  (SZ,SV).  Since  this  system  is  overconstrained,  we  solve  it  using 
the  least-squares  method  by  solving  the  2x2  system 


ArAA  =  ~ATb 


(6) 


where  A  is  the  2x3  matrix,  A  =  (£i,6y),  and  6  the  righthand  side  vector. 

If  the  singularity  has  higher  order,  the  system  (6)  is  singular.  In  this  case 
we  determine  which  higher  order  partials  also  vanish.  For  each  vanishing 
higher  order  partial  derivative  h  of  /,  the  matrix  A  is  augmented  by  the  row 
{hz(Pi),hv(Pi))  and  the  vector  6  by  the  entry  h(Pj).  This  process  continues 
until  D  =  At A  has  full  rank.  For  exceptional  values  it  is  possible  that  D 
has  rank  1  or  zero  even  though  no  additional  partials  of  /  vanish.  This 
means  that  we  happen  to  approach  the  singular  point  crossing  a  specific 
algebraic  curve  given  by  the  symbolic  determinant  of  D.  In  this  case,  a 
random  perturbation  of  the  point  should  correct  the  problem.  So  far,  we 
have  not  encountered  this  problem  in  practice. 


5.3  Passing  to  the  Transformed  Curve 

By  a  translation  of  /,  the  coordinate  system  is  centered  at  the  singular 
point  q  just  found.  This  may  introduce  spurious  terms  that  are  controlled 
based  on  the  information  obtained  during  the  iteration  locating  q.  Recall 
that  the  vector  b  in  the  iteration  contains  all  partial  derivatives  that  vanish. 
Consequently,  if  the  partial  h  appears  in  b,  then  the  corresponding  monomial 
term  must  be  absent  in  the  translated  curve.  For  example,  let  /  be  the 
translation  of  /,  and  assume  that  b  contains  the  vanishing  partials  /*,  /„, 
/**>  fvv >  /* y>  /zyy>  and  /yyy  Then  /  must  not  contain  the  terms  i,  y,  x2, 
y2,  xy,  xy2,  and  y3.  Should  such  terms  appear  in  /  with  small  coefficients, 
due  to  numerical  imprecision,  they  are  now  removed. 

Having  centered  the  singularity  at  the  origin,  we  apply  the  quadratic 
transformation  (5).  Since  this  transformation  maps  the  line  x  =  0  to  infinity, 
the  branch  of  /  we  traverse  must  not  have  the  y-axis  as  tangent.  If  it  does, 
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we  rotate  /  by 


x'  —  y  y'  =  -x 

before  applying  the  quadratic  transformation. 

In  practice,  one  applies  instead  the  quadratic  transformation 

x 

xi  =  -  yi  =  y 

y 

in  which  the  i-axis  becomes  the  exceptional  line.  Which  quadratic  trans¬ 
formation  is  used  is  decided  based  on  the  current  tangent  direction  of  /  in 
the  traversal. 

5.4  Branch  Orientation 

We  give  a  standard  orientation  to  the  curve  /( x,y)  =  0  by  the  tangent  vector 
(“ /v>/z)-  This  orientation  is  not  intrinsic  in  the  sense  that  -f{x,y)  =  0 
is  the  same  curve  but  with  opposite  standard  direction.  Given  a  consistent 
traversal  direction  of  a  branch,  we  observe  that  the  standard  direction  of 
the  curve  may  reverse  at  certain  singularities.  For  example,  the  orientation 
of  /  =  y2  -  x2  -  i3  =  0  is  as  shown  in  Figure  5.3.  Consequently,  when 
traversing  the  curve  from  p  to  q,  we  first  move  in  the  standard  direction, 
but  after  the  singularity  we  move  in  the  opposite  direction.  Figure  5.4  shows 
that  this  reversal  does  not  happen  at  all  singularities. 

Geometrically,  the  apparent  orientation  reversal  is  understood  when  con¬ 
sidering  f(x,y)  =  0  as  the  intersection  curve  of  /(*, y)  -  z  —  0  and  z  =  0. 
Along  the  intersection  curve,  the  projection  onto  the  x-y  plane  of  the  sur¬ 
face  gradient  -1)  is  just  the  curve  normal.  Thus  the  normal  reversal 

that  causes  the  changed  standard  orientation  is  merely  a  rotation  of  the 
.urface  normal  in  3-space.  Whether  a  branch  suffers  this  reversal  depends 
on  the  global  topology  of  the  singularity.  Briefly,  the  orientation  reverses 
if  the  branch  intersects  an  odd  number  of  other  branches,  with  the  proper 
definition  of  intersection.  Nevertheless,  a  local  correspondence  between  the 
orientation  of  /  and  its  proper  transform  g  can  be  established  outside  the 
singularity  from  which  we  can  deduce  whether  the  standard  orientation  has 
reversed,  without  having  to  analyze  the  topology  of  the  singularity. 

Let  p  =  (ao,6o)  be  a  nonsingular  point  of  /,  where  ao  #0.  To  p  corre¬ 
sponds  the  point  pi  =  (ao,i>o/oo)  of  the  transformed  curve  g.  Centered  at 
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p,  the  curve  /  has  the  place 


x(s )  =  UO  +  als  +  02S1  +  •  •  • 
y{s)  —  bo  +  bis  +  &2fi2  +  •  •  • 

and  centered  at  pi ,  the  curve  g  has  the  place 

X!(s)  =  x(s)  =  a0  +  ais  +  a2S2H - 

S/l(s)  =  Co  +  CjS  +  C2S2  +  •  ■  • 

We  assume  that  the  traversal  of  /  at  p  proceeds  by  increasing  value  of  s.  Note 
that  the  traversal  direction  need  not  agree  with  the  standard  orientation 
(—/„,/*).  Since  x(s)  =  xi(s),  the  curve  and  its  transform  are  oriented  the 
same  way,  and  traversing  g  by  increasing  s  is  equivalent  tc  traversing  /  by 
increasing  s. 

Since  t/x(s)  =  y(s)/x(s),  we  divide  the  two  power  series  and  compare  the 
resulting  coefficients  with  the  c*.  We  obtain 

co  =  bQ/ ao 

Mo  ~  Mo 

Cl  =  - ? -  (7) 

°o 


Since  /  is  not  singular  at  p,  g  is  not  singular  at  p\.  Hence  both  curves  have 
a  Taylor  series  at  these  points  so  that  ai  is  proportional  to  —fv  and  to  —  gy, 
while  61  is  proportional  to  fx,  and  ci  is  proportional  to  gx.  We  now  obtain 
from  (7)  that 


9v 

9z 


<*fv 

xfx  +  y  fy 


So,  if  we  relate  the  traversal  direction  of  /  to  the  standard  orientation 
(~fv,fx),  then  a  relates  the  corresponding  traversal  direction  of  g  to  the 
standard  orientation  {-gv,gx)  of  g,  and  vice  versa.  Since  the  fully  desingu- 
larized  branch  cannot  experience  an  orientation  reversal,  we  have  a  method 
to  maintain  consistent  traversal  direction  through  singularities. 
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5.5  Implementation 


The  planar  curve  traversal  has  been  implemented  in  Lisp  on  a  Symbolics 
Lisp  machine.  A  prototype  was  previously  implemented  by  C.  Bajaj  on  a 
VAX  8600.  Figure  5.4  shows  an  example  of  a  traversal  requiring  iterated 
desingularization,  as  well  as  the  traversals  along  the  proper  transforms  in 
the  vicinity  of  the  singularity.  For  simplicity,  the  singularity  was  already 
positioned  at  the  origin. 


6  Singularities  on  Space  Curves 

Every  algebraic  space  curve  is  the  birational  image  of  an  algebraic  plane 
curve.  It  follows  that  the  singularities  of  space  curves  are  not  qualita¬ 
tively  different  from  those  of  plane  curves.  Two  possibilities  exist  for  tracing 
through  space  curve  singularities: 

1.  Construct  a  birational  map  from  the  given  space  curve  to  a  planar 
curve,  trace  the  planar  curve,  and  lift  the  resulting  points. 

2.  Desingularize  the  space  curve  directly. 

Since  the  intersection  curve  in  general  has  degree  equal  to  the  product 
of  the  surface  degrees,  the  birationally  equivalent  plane  curve  must  have 
high  degree  and  may  be  computationally  less  tractable.  Working  with  the 
space  curve  directly  is  therefore  more  attractive.  However,  desingularizing 
the  space  curve  directly  must  address  the  fact  that  the  curve  is  given  as  the 
intersection  of  two  surfaces.  If  we  are  to  work  with  this  representation,  then 
we  need  to  apply  quadratic  transformations  that  desingularize  the  curve 
and  the  intersecting  surfaces  as  well.  Hence  the  approach  to  space  curve 
desingularization  found  in  the  standard  literature,  e.g.  [5],  will  not  work,  and 
more  research  is  required  to  work  out  techniques  suitable  for  computation. 
We  restrict  our  discussion  therefore  to  the  method  of  reducing  the  space 
curve  to  a  plane  curve. 

The  simplest  way  to  map  the  space  curve  to  a  planar  curve  is  by  projec¬ 
tion.  Orthographic  projection  along  a  principal  axis  is  done  by  elimination 
of  a  variable,  using  resultants.  Other  directions  require  a  rotation  of  the 
coordinate  system  prior  to  projection.  For  space  curve  singularities  where 
at  least  one  of  the  surfaces  has  a  nonzero  gradient,  orthographic  projec¬ 
tion  onto  the  tangent  plane  is  conceptually  ideal.  The  major  computational 
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problem  would  be  the  inefficiency  of  the  resultant  computation  for  surfaces 
of  high  degree.  Moreover,  branches  intersecting  the  surface  normal  above 
or  below  the  tangent  plane  are  also  projected  and  unnecessarily  increase  the 
complexity  of  the  singularity. 

A  different  method  to  map  a  space  curve  to  the  plane  is  to  find  a  rational 
surface  containing  the  curve,  parameterizing  this  surface,  and  then  substi¬ 
tuting  the  parametric  equations  into  one  of  the  implicit  surface  equations, 
say  g.  This  method  has  the  advantage  of  treating  all  singularities,  not  only 
those  at  a  nonzero  surface  gradient.  We  describe  the  method  below  and  give 
several  examples.  It  has  not  been  implemented  yet. 

6.1  Monoid  and  Cone  Representation 

It  is  well  known  that  every  algebraic  space  curve  /  n  g  can  be  represented  as 
the  intersection  of  a  monoid  and  a  cone,  [13].  A  monoid  is  a  rational  surface 
of  degree  m  that  contains  an  m  —  1  fold  point.  Simple  examples  include 
all  planes,  quadrics,  and  the  Steiner  surface.  When  the  m  -  1  fold  point  is 
brought  to  the  origin,  the  monoid  equation  takes  the  form 

wHm-i(x,y,z)  +  Hm(x,y,z)  =  0 

where  /fm_j  is  homogeneous  of  degree  m  -  1  and  Hm  is  homogeneous  of 
degree  m. 

We  will  be  interested  in  determining  the  monoid  containing  a  given  space 
curve  /  n  g  and  its  parametric  representation.  We  will  not  determine  the 
cone,  since  it  is  not  needed.  Moreover,  the  parameterization  of  the  monoid 
is  incompatible  with  the  cone  equation.  The  procedure  for  determining  the 
monoid  is  based  on  the  projective  form  of  the  surfaces  /  and  g.  As  in  [12], 
we  proceed  as  follows. 

First,  homogenize  f(x,y,z)  and  g{x,y,z)  so  as  to  obtain  F{w,x,y,z) 
and  G(w,x,y,  z).  As  long  as  u;  ^  0,  the  curve  F  n  G  is  identical  to  /  n  g. 
We  select  one  of  the  base  points  of  the  projective  coordinate  system,  say 
(1,0, 0,0),  as  the  m  -  1  fold  monoid  point.  This  implies  using  u;  as  the 
main  variable  in  the  computation  below.  The  base  point  (0, 1,0,0)  would 
correspond  to  selecting  x  as  the  main  variable,  and  so  on. 

We  write  both  F  and  G  as  polynomials  in  the  main  variable,  w, 

=  Unwn  +  +  ••  •  +  UlW  +  Uo 

=  vniwn>  4-  Vui-itu”'-1  -I - hvitx/  +  vo 


F 

G 
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Without  loss  of  generality  we  assume  that  n  >  n'  >  1.  We  compute  the 
polynomials 

F\  =  unu)n~n>G  -  vn>F 
G\  =  {uqG  —  v0F)/w 

Note  that  both  F\  and  G i  contain  the  intersection  curve  of  F  and  G. 

Both  F\  and  G i  have  degree  at  most  n  -  1  in  to.  If  one  of  them  is 
linear  in  to,  then  we  stop;  we  have  found  the  monoid  equation.  If  neither  is 
linear,  then  we  repeat  the  calculation  using  Fi  and  Gi  in  place  of  F  and  G. 
Since  at  each  step  the  maximum  degree  in  to  is  lowered  at  least  by  one,  the 
computation  derives  the  monoid  equation  after  at  most  n  steps  in  the  form 

wHm-i(x,y,z)  +  Hm(x,y,z)  =  0 

The  monoid  is  parameterized  by  intersecting  it  with  lines  through  the  m  —  1 
fold  point.  Let  a  :  b  :  c  be  the  direction  ratios  of  these  lines,  then  the  monoid 
is  parameterized  by 

w(a,  b,  c)  =  -Hm-i(a,b,c)/Hm{a,b,c) 
x(a,b,c)  =  a 
y(a,b,c)  =  b 
z(a,b,c )  =  c 

This  parameterization  is  projective,  that  is,  (o,6,c)  are  the  coordinates  of  a 
two-dimensional  projective  parameter  space. 

The  parametric  forms  are  now  substituted  into  the  equation  of  G  and 
give  the  desired  plane  curve,  in  homogeneous  form. 

6.2  Examples 

We  illustrate  the  method  with  two  examples.  First,  consider  the  intersection 
curve  of  the  cylinder  F  =  x2  +  z2  +  2 zw  =  0  and  the  sphere  G  =  z2  +  y2  + 
z 2  +  4 zw  =  0.  The  intersection  curve  is  an  irreducible  degree  4  space  curve 
with  a  nodal  singularity  at  the  origin  shown  in  Figure  6.1. 

The  cylinder  is  a  monoid  with  the  m  - 1  fold  point  at  the  origin  (1, 0, 0, 0). 
Since  the  point  of  interest  on  the  space  curve  is  the  origin,  we  determine 
a  different  monoid  whose  m  -  1  fold  point  is  not  the  origin.  We  choose 
(0,0,0, 1),  making  z  the  main  variable.  Accordingly,  we  compute 

Fi  =  G  -  F  =  y2  +  2zw 
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The  parameterization  of  Fi  is  then 


z  =  - 2a/c 2 

w  =  a 

x  =  b 

y  =  c 

Substitution  into  G  yields  the  plane  curve 

64  +  4a2  (c2  -  62)  =  0 

Dehomogenizing  with  a  —  1  yields  64  —  4 (e2  -  62)  =  0.  This  curve  is  shown 
in  Figure  6.2. 

As  a  more  complicated  example,  consider  the  intersection  of  the  torus 
(x2  +  y2  +  z2  -  u/2)2  +  8u;2 (z2  -  x2  -  y2  -  u>2)  +  16ti;4  =  0  with  the  ellipsoid 
36(z  -  tu)2  +  4(y  -  w)2  +  9x2  -  36u>2  =  0.  The  monoid  computation,  as 
described,  yields  a  surface  of  degree  12  in  4  steps.  Its  equation  contains  an 
extraneous  factor  of  degree  4.  Substitution  into  the  ellipsoid  equation  thus 
yields  a  plane  curve  of  degree  16  that  factors  into  a  degree  2  component,  a 
degree  6  component,  and  a  degree  8  component. 

A  better  solution  to  the  problem  is  to  parameterize  the  ellipsoid  since 
it  is  a  monoid.  To  do  so,  we  first  translate  the  coordinate  system  to  the 
point  (1,0, 1,2),  which  is  on  the  ellipsoid,  and  parameterize  with  w  as  main 
variable.  This  results  in  a  plane  curve  of  degree  12  that  could  not  be  factored 
by  Macsyma. 

6.3  Remarks 

Only  in  the  plane  curve  case  do  we  know  of  a  simple  algebraic  procedure 
achieving  complete  desingularization  at  minimal  computational  cost.  The 
projection  of  space  curves  to  planar  curves  seems  to  require  significant  ma¬ 
chinery,  for  determining  the  monoid  equation,  and  for  eliminating  extrane¬ 
ous  components  introduced  in  the  process.  In  the  case  of  monoid/cone  in¬ 
tersection  these  extraneous  components  are  all  lines,  hence  would  be  simple 
to  exclude.  The  cone  is  determined  by  a  resultant  computation,  and  yields 
a  homogeneous  polynomial  in  three  variables.  Unfortunately,  its  equation 
is  unsuitable  for  substituting  the  monoid  parameterization.  In  simple  cases, 
there  exist  certain  reparameterizations  that  circumvent  this  problem. 
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If  one  of  the  surfaces  /  or  g  is  known  to  be  rational  it  can  be  advanta¬ 
geous  to  parameterize  it  directly.  In  the  case  of  quadrics  this  amounts  to  a 
coordinate  translation  that  brings  one  of  the  base  points  of  the  coordinate 
system  onto  the  surface.  In  the  case  of  cubics,  [2]  gives  parameterization 
algorithms.  Many  other  surfaces,  including  the  torus,  are  also  rational  with 
known  standard  parameterizations.  Direct  parameterization  side-steps  a  po¬ 
tentially  lengthy  monoid  derivation.  It  should  be  noted,  however,  that  the 
resulting  plane  curve  is  not  necessarily  of  minimum  degree. 

It  appears  that  there  are  simple  quadratic  transformations  of  space  that 
achieve  surface/surface  intersection  desingularization  much  as  in  the  planar 
case.  More  research  is  needed  to  explore  the  exact  effect  of  these  transfor¬ 
mations,  and  the  issues  involved  in  realizing  them  by  computation. 
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Figure  5.5 


Recursive  Desingularization 
2x4  -  3x2y  +  y2-  2ys  +  y4 


2i2  -  3y!ix  +  y\  -  2yfxx  +  y\x\ 


2  -  3y2  +  y\  -  2y\x\  +  y\x\ 


Figure  6.1 


Cylinder  -  Sphere  Intersection 
x2  +  y2  +  z1  +  4z  n  x2  +  z1  +  2  z 


Figure  6.2 

Corresponding  Plane  Curve 
b*  -  4(c2  -  b 2) 


